package edu.gmu.atelier;

/**
 * Assuming A is an orthogonal matrix, computes b* and x* using 
 * the normal projection equation bHat = A x*.
 * <p>
 * The normal equation A^T A x = A^T b, I believe assuming A is orthogonal
 * then this becomes x = A^T b.
 * 
 * I NEED TO LOOK AT THIS MORE - THIS IS SPECIAL CASE, NEED GENERAL CASE
 * SO LOOK AT SOLVING NORMAL EQUATION AND xHat = R^-1 Q^T.
 * 
 * @author James H. Pope
 */
public class LSOrthogonalSolver implements Solver
{
    private Matrix a    = null;
    
    public LSOrthogonalSolver( Matrix a )
    {
        this.a = a;
        
        // Must be orthogonal
        validate(a);
    }
    
    //------------------------------------------------------------------------//
    // Interface methods
    //------------------------------------------------------------------------//
    public Matrix getA()
    {
        return this.a;
    }
    
    public Vector solve( Vector b )
    {
        // Note: Orthogonal A^T A = D, orthonormal A^T A = I
        int n = a.getCols();
        int m = a.getRows();
        Vector xStar = a.getFactory().makeVector(n);
        for (int i = 0; i < n; i++)
        {
            Vector qi    = a.getVector(i);
            Element num  = b.dot(qi);
            Element den  = qi.dot(qi);
            Element coeff= num.mult( den.inverse() ); // Note: Weight in xStar
            xStar.setElement(i, coeff);
            //Vector projV = vi.dot( coeff );
            //bStar = bStar.add( projV ); // overly exspensive imho
        }
        return xStar;
    }
    
//    private void decompose()
//    {
//        // Must be orthogonal
//        validate(a);
//        
//        int n = a.getCols();
//        int m = a.getRows();
//        this.bStar = a.getFactory().makeVector(m);
//        this.xStar = a.getFactory().makeVector(n);
//        for (int i = 0; i < n; i++)
//        {
//            Vector vi    = a.getVector(i);
//            Element num  = b.dot(vi);
//            Element den  = vi.dot(vi);
//            Element coeff= num.dot( den.inverse() ); // Note: Weight in xStar
//            this.xStar.setElement(i, coeff);
//            Vector projV = vi.dot( coeff );
//            bStar = bStar.add( projV ); // overly exspensive imho
//        }
//    }
    
    private static void validate(Matrix a)
    {
        if( a.isOrthogonal() == false)
        {
            String s = "Matrix not orthogonal, cannot perform least squares";
            throw new IllegalArgumentException(s);
        }
    }
    
    
    public static Vector computeXstar( Matrix a, Vector b )
    {
        // Must be orthogonal
        validate(a);
        
        int n = a.getCols();
        int m = a.getRows();
        Vector xHat = a.getFactory().makeVector(n);
        for (int i = 0; i < n; i++)
        {
            Vector vi    = a.getVector(i);
            Element num  = b.dot(vi);
            Element den  = vi.dot(vi);
            Element coeff= num.mult( den.inverse() ); // Note: Weight in xHat
            xHat.setElement(i, coeff);
            //Vector projV = vi.dot( coeff );
            //bHat = bHat.add( projV ); // overly exspensive imho
        }
        return xHat;
    }
    
    public static Vector computeBstar( Matrix a, Vector b )
    {
        // Must be orthogonal
        validate(a);
        
        int n = a.getCols();
        int m = a.getRows();
        Vector bHat = a.getFactory().makeVector(m);
        for (int i = 0; i < n; i++)
        {
            Vector vi    = a.getVector(i);
            Element num  = b.dot(vi);
            Element den  = vi.dot(vi);
            Element coeff= num.mult( den.inverse() ); // Note: Weight in xHat
            Vector projV = vi.mult( coeff );
            bHat = bHat.add( projV ); // overly exspensive imho
        }
        return bHat;
    }
    
}
